The interplay between boundary conditions and flow geometries in shear banding: 
hysteresis, band configurations, and surface transitions 
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We study shear banding flows in models of wormlike micelles or polymer solutions, and explore 
the effects of different boundary conditions for the viscoelastic stress. These are needed because the 
equations of motion are inherently non-local and include "diffusive" or square-gradient terms. Using 
the diffusive Johnson-Segalman model and a variant of the Rolie-Poly model for entangled micelles 
or polymer solutions, we study the interplay between different boundary conditions and the intrinsic 
stress gradient imposed by the flow geometry. We consider prescribed gradient (Neumann) or value 
(Dirichlet) of the viscoelastic stress tensor at the boundary, as well as mixed boundary conditions 
in which an anchoring strength competes with the gradient contribution to the stress dynamics. We 
find that hysteresis during shear rate sweeps is suppressed if the boundary conditions favor the state 
that is induced by the sweep. For example, if the boundaries favor the high shear rate phase then 
hysteresis is suppressed at the low shear rate edges of the stress plateau. If the boundaries favor the 
low shear rate state, then the high shear rate band can lie in the center of the flow cell, leading to 
a three-band configuration. Sufficiently strong stress gradients due to curved flow geometries, such 
as that of cylindrical Couette flow, can convert this to a two-band state by forcing the high shear 
rate phase against the wall of higher stress, and can suppress the hysteresis loop observed during a 
shear rate sweep. 



I. INTRODUCTION 

Microscopic models for certain viscoelastic fluids, such 
as wormlike micelles [l[ and high molecular weight poly- 
meric liquids [HEl predict unstable homogeneous station- 
ary states for controlled average shear rates. Specifically, 
the coupled equations of motion for the fluid flow and 
microstructural quantities, such as the viscoelastic stress 
£ or orientational order, predict a nonmonotonic con- 
stitutive curve, i.e. the steady state total shear stress 
T xy (or equivalently the applied torque) as a function of 
applied shear rate 7 for homogeneous flows. This con- 
stitutive curve typically is multivalued with a region of 
decreasing shear stress as a function of increasing shear 
rate, which is a hallmark of hydrodynamic instability [3] . 
Understanding this nonlinear behavior is important for 
practical applications such as injection moulding of plas- 
tics and drilling muds used in bore holes Q ■ The simplest 
resolution of the instability is an inhomogeneous state 
with macroscopic regions of high and low shear rates, 
known as shear bands [f| 0L which has been widely ob- 
served in wormlike micelles [1, [t| Qi| fill . lamellar surfac- 
tant solutions [l2j j and liquid crystals jl3l| . The resulting 
experimental signature of this inhomogeneous state is a 
stress plateau as a function of controlled average shear 
rate. We denote this experimentally determined relation 
between total shear stress and average shear rate as the 
flow curve; in shear banding flow this curve incorporates 
inhomogeneous flows and is distinct from the constitutive 



curve, which cannot be measured because of the consti- 
tutive instability. In Fig.Q]the dashed line shows the con- 
stitutive curve while the triangles show flow curves that 
could be measured upon either increasing or decreasing 
average shear rate ramps. 

Experimentally, it is clear that the geometry of the 
shear cell affects the positioning of the shear bands. The 
simplest structure is seen in the cylindrical Couette ge- 
ometry where two bands form, with the higher shear rate 
band near the inner wall (rotor) [l4l |. In the cone and 
plate geometry a different band configuration has been 
reported: two low shear rate regions next to the cone 
and plate separated by a high shear rate region [§]. One 
explanation offered for that result was that secondary 
flows stabilized the center band [l5j]. We will show that 
a boundary condition that prescribes a value for the poly- 
meric stress tensor similar to that of the low shear rate 
phase can induce this three-band configuration. 

The experimental rheological features of shear band- 
ing have been studied in detail for wormlike micellar so- 
lutions [l(| El, [HI- A constitutive curve such as the 
dashed line in Fig. Q] can only be inferred because the 
negative slope region is unstable; however, the exper- 
imentally measured flow curve, indicated by triangles, 
typically shows a stress plateau which spans the range of 
average applied shear rates at which the system shows 
shear banding. During shear rate ramps from rest, hys- 
teresis is often observed at the start of this plateau: the 
stress increases past the plateau stress and follows the 
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constitutive curve until, at a time that depends on the 
rate of the ramp, a high shear rate band develops and 
the total stress decreases onto the steady state plateau 
[ToL [Til. []~7l] . However, if the shear rate is reduced from a 
point on the plateau, the low shear rate branch of the con- 
stitutive curve is intersected directly by the flat plateau. 
This behavior is reminiscent of nuclcation and hysteresis 
at first order phase transitions, and we will show below 
that the boundary conditions can influence this hystere- 
sis, in a manner analogous to heterogeneous nucleation. 
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FIG. 1: Dashed line: the constitutive curve (dimensionless 
specific torque T as a function of dimensionless shear rate 7) 
for the diffusive Johnson- Segalman (DJS) model in cylindri- 
cal Couette flow. The triangles show simulated steady state 
flow curves for zero gradient (Neumann) boundary conditions, 
for increasing (>) and decreasing (<) the average shear rate 
(7) from small or large average shear rates respectively. The 
model parameters are a — 0.3, e = 0.05, T> = 10~ 5 , and the 
gap size is 0.1% of the radius (q — 0.005). The physical differ- 
ence between the two bands is illustrated by ellipsoids whose 
principal axes' directions and length ratios coincide with those 
of the viscoelastic stress tensor S_l and YIh in, respectively, 
the low and high shear rate phases. 

In modelling shear banding systems such as wormlikc 
micelles the most important degrees of freedom are the 
fluid velocity v and a viscoelastic contribution to the 
stress, S; here we neglect other other potentially impor- 
tant quantities such as concentration and micellar length. 
Among models with non-monotonic constitutive curves 
for homogeneous flow, spatially local ones such as the 
Johnson-Segalman model in its original form [l8j show 
high sensitivity to the flow history and do not give the 
well-defined and unique stress plateau observed in experi- 
ments |19| . In contrast, if a non-local diffusive term (typ- 
ically proportional to V 2 S) is included in the constitu- 
tive equation for the viscoelastic stress (e.g. the diffusive 
Johnson-Segalman (DJS) model), there is only a single 
value of the total stress for which the interface between 



bands is stable and stationary [20[. The diffusive term 
arises from microscopic physical mechanisms such as the 
diffusion of molecules that carry stress [2l|, the persis- 
tence length of wormlike micelles (22j . or hydrodynam- 
ics [23L |24| . The resulting equation governing inhomoge- 
neous steady states is thus a spatial differential equation 
for the viscoelastic stress, which necessitates a bound- 
ary condition that is, at present, unknown. The most 
frequently used boundary condition on the viscoelastic 
stress S has been that of zero stress gradient parallel to 
the boundary normal [20| . This boundary condition pre- 
dicts the flow curve shown in Fig. [1] for the DJS model. 
As noted above, this flow curve shows hysteretic behav- 
ior at the start of the stress plateau similar to that seen 
in experiment [TTI ] . with a stress overshoot during an in- 
creasing shear rate ramp. 

Other theoretical studies of shear banding have used a 
fixed value of the polymeric stress (2f| (2(| [23 at the wall 
(Dirichlet boundary conditions). Cook and Rossi con- 
sidered a two-fluid model for wormlike micelles in which 
the micelles were assumed to align at the wall parallel to 
the flow direction, in both planar and cylindrical Cou- 
ette flow. They found that the flow curve deviates sig- 
nificantly from the constitutive curve near the low shear 
rate branch (2f| [26| , exhibited hysteresis only at the high 
shear rate side of the stress plateau, and possessed a 
three-band state with the boundary condition induced 
high shear rate bands near the wall. Qualitatively sim- 
ilar results were found by Picard and co-workers in a 
scalar model for shear banding in a yield stress solid [13] • 

In this paper we perform a more detailed study of the 
effects of different boundary conditions for the viscoelas- 
tic stress at the wall, So- We consider strong "anchoring" 
in which Dirichlet conditions apply, and interpolate be- 
tween oblate and prolate forms of So- We also study 
mixed boundary conditions, in which So is determined 
by a balance between spatial gradients and surface an- 
choring. Moreover, we address the important interplay 
between the boundary conditions and intrinsic inhomo- 
geneity of the flow determined by the flow geometry (to 
compare cylindrical and planar Couette flow with cone 
and plate flow, for example). 

The organization of the paper is as follows. In Sec- 
tion [n] we outline the details for solving the equations of 
motion within the creeping flow approximation in cylin- 
drical Couette flow, and describe two constitutive equa- 
tions for the viscoelastic stress: the DJS model [2(| and a 
tube-based model which has been developed for polymer 
solutions and wormlike micelles [ID, H^]. The physical 
motivation for the boundary conditions is discussed in 
Section IIIIl After reviewing previous results for Neu- 
mann boundary conditions in Section [Vj the main new 
results are presented in Section [VTl We organize these as 
follows: (1) the effect of different fixed value (Dirichlet) 
boundary conditions on the flow curves and hysteresis 
upon increasing and then decreasing the average shear 
rate; (2) the interplay between boundary conditions and 
the stress gradient imposed by the flow geometry (e.g. 
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cylindrical Couette flow) and their effects on the stable 
position of the shear bands; and (3) the effect of mixed 
boundary conditions, in which the viscoelastic stress at 
the surface is influenced by both the bulk constitutive re- 
lation and the flow geometry. We demonstrate the possi- 
bility of a transition between different effective boundary 
conditions as a function of the anchoring strength. 

II. EQUATIONS OF MOTION 

A. Creeping Flow Approximation 

In a wormlike micellar system we assume that the total 
stress T can be separated into contributions from the 
Newtonian solvent and a viscoelastic stress 53 from the 
micelles: 

T = 217D + 53 - pi. (11.1) 

Here I denotes the identity matrix, ?/ is the solvent shear 
viscosity, D = \ [Vv + (Vv) T ] is the symmetric veloc- 
ity gradient tensor, and p is the isotropic pressure deter- 
mined by incompressibility (V-v = 0). We are interested 
here in the creeping flow regime (low Reynolds number), 
for which the force balance has the form 

V • T = 0. (II.2) 

We perform calculations explicitly for Couette flow be- 
tween concentric cylinders and assume unidirectional 
flow v = v(r,t)0, where r and 6 are the usual cylindrical 
coordinates. In this case, Eqs. pi. 11111. 2j) imply that 

rrr(r,t) + Z r0 (r,t) = ^, (II.3) 

where 

7(r,t)=4^, (0.4) 
or r 

and r, the specific torque per unit cylinder height per 
radian, is a constant of integration. We assume no-slip 
boundary conditions on the velocity field. For cylindrical 
Couette flow with a fixed outer cylinder and a rotating 
inner cylinder of velocity V, the no-slip boundary condi- 
tions lead to the following global constraint 

V f Rl dr 
Ri Jr 2 r 
which implies that the specific torque T is given by 

r flj-fli _JnV f R2 y dr\ 
T -RM- 2 {lh-J Rl ^-)- (IL6) 

To simplify the equations we define the spatial variable 

old 



where q = In |p. Making use of this variable change, the 
creeping flow equation becomes 

r?7 = Te- 2 ^ - S3 r0 , (II.8) 

with the specific torque is given by 

r = T ^^((^)-'7(7)) I (0.9) 

where the spatial averages are simplified to ((•)) = 
/ (-)dx. This expression can be used to calculate T for an 
imposed average shear rate, and hence render Eq. (|II.8p 
an integral equation determining the relation between the 
local and average shear rates and viscoelastic shear stress. 
A complete description of the dynamics requires an equa- 
tion of motion for the viscoelastic stress 53, examples of 
which are described below. 



B. Constitutive Equations 

1. Diffusive Johnson-Segalman (DJS) Model 

One of the simplest tensorial models that produces 
a non-monotonic constitutive curve is the DJS model, 
which has been studied in detail [U, ES HI HI ■ Phys- 
ically it is motivated by neglecting all but the lowest 
modes of vibration of polymer chains, thereby represent- 
ing them as elastic dumbbells with span R [3l|. The 
viscoelastic stress tensor 53 is related to the extension of 
the dumbbells by 53 = G{— I + fe(RR)), where expres- 
sions for the normalization k and the modulus G can be 
derived in terms of the parameters in the underlying mi- 
croscopic polymeric models. Here ((•)) denotes a thermal 
average. The DJS model [2(| [HJ provides the constitu- 
tive equation for the evolution of the viscoelastic stress: 

S +-53 = 2-D + PV 2 53, (11.10) 

T T 

where 

£= (<9 t +v-V)53+(0-53-53-0)-a(D-53+53-D), (11.11) 

t is a relaxation time, the "polymer" viscosity (i is related 
to the modulus by G = (i/t, and fi = | [Vv — (Vv) T ] . 
The total stress comprises the viscoelastic stress of the 
DJS model and a Newtonian contribution, according to 
Eq. (|II.1[) . and the viscosity ratio e = rj/fj, controls the 
balance between the two stresses. The 'slip parameter' 
a, which describes non-affine stretch of the dumbcll with 
respect to the extension of the flow, allows for a non- 
monotonic constitutive curve for < \a\ < 1 and e < 1/8. 

The "diffusion" term 2?V 2 53 describes non-local relax- 
ation of the viscoelastic stress and is necessary to describe 
strongly inhomogeneous flow profiles [2(| • Because of this 
term the steady shear banding state obeys a spatial dif- 
ferential equation, which must solved subject to bound- 
ary conditions specified at the walls of the flow cell. The 
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solvability condition for a stationary interface leads to 
a unique total shear stress plateau for imposed average 
shear rates in the non-monotonic portion of the constitu- 
tive curve [24| • The characteristic width I of the interface 
between shear bands is given by I — \JT>t . 

For convenience we define the modulus G = /u/r. In 
cylindrical Couette flow the constitutive equation has the 
following components [20| : 



C T, r e 
where 



2T>e~ 2qx 
-(1 - a)jE re + — j-2— (S fl 

(1 + o)7£ re 2 (S fl 



G 



1 -a. 



1 + a, 



-S rr ) (II.12a) 
-S rr ) (II.12b) 

(II.12c) 



9 1 Pe" 2 ^ d 2 

£ = 1 , 

dt t q 2 R 2 dx 2 



(11.13) 



2. Reptation- Reaction and Rolie-Poly Models 

In Cates' model [33| of wormlike micellar solutions the 
micelles are assumed to relax in a tube, and the devia- 
toric viscoelastic stress is defined by £ = G{— I+3(uu}), 
where the unit vector u describes the local orientation 
of tube segments. The viscoelastic stress is given by a 
history integral over the second moment of the correla- 
tion function (uu) , weighted by the distribution of tube 
segments as they are created and subsequently vacated 
by the wormlike micelles. The resulting integral equation 
can be re- written as the following differential equation by 
using a decoupling approximation to remove fourth order 
moments 12911: 



v 
£ 



1. 



2GD--E : Vv (1+ (1 



/3)[£/G])+2?V 2 £, 
(11.14) 



where 

v 
£= 



(9 t +v-V)£+(n-£-£-f2)-(D-£+£-D). (11.15) 



and (3 = 0. A stress diffusion term, which did not appear 
in the original formulation, has been included here. The 
non-linear term in this equation preserves the traceless 
property of the deviatoric stress. 

Interestingly, for (3 ^ this model is identical to an 
extension of the original Doi-Edwards (DE) theory of en- 
tangled polymers |2J , which has a constitutive instability 
and therefore predicts shear banding. This extension to 



the DE theory incorporates the enhanced release of poly- 
mer entanglements due to convection (3~ij|, which increases 
the polymer stress and can "cure" the DE instability for 
sufficiently strong convected constraint release (CCR). 
CCR and tube stretching were incorporated in Q to de- 
scribe polymer melts and wormlike micelles, and later 
simplified to a differential version called the Rolie-Poly 
model [23| . In the version of the Rolie-Poly model used 
in Eq. (|II.14[) the tube length is assumed to be relaxed 
("non-stretching"). The CCR parameter (3 (0 < (3 < 1) 
corresponds to the frequency of the release of polymer 
entanglement constraints due to convection by the flow, 
and Ref. [H| used (3 = 1 to model a well-entangled poly- 
mer melt without a constitutive instability. For small 
enough (3 a constitutive instability results. 
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FIG. 2: Regions of parameter space f3 and e for which 
the DRP model has either a monotonic (unshaded) or non- 
monotonic (shaded) constitutive curve. 



The DE constitutive instability is thought to be an un- 
physical prediction of polymer melts [|[ , but a series of re- 
cent experiments have reignited the interest in this insta- 
bility by showing data, such as velocity profiles, consis- 
tent with shear banding in polymer solutions [35], [36l. |37|] . 
Hence, the Rolie-Poly model is an alternative simple con- 
stitutive model to the DJS model for wormlike micelles, 
and by tuning the parameter (3 it can be used to address 
polymer solutions and melts. 

As with the DJS model, we define the modulus by 
G = [i It and the ratio of Newtonian (or "solvent") to 
viscoelastic stress by e = rj/fj,. For a given e a non- 
monotonic constitutive curve is found if the degree of 
convected constraint release (3 is small enough (Fig. [2]). 
We will refer to the model described above as the diffusive 
Rolie-Poly (DRP) model. In cylindrical Couette flow the 
components of the DRP model evolve according to 
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2 2Ve~ 2qx 

C £ rr = --7£ re [1 + (1 + /3)(S rr /G)] + 2 - S rr ) (II.16a) 

2 2Ve~ 2qx 

C S ee = --7£ r <? [1 + (1 + P)(?ee/G)] + {^ee - Z rr ) (II.16b) 

2 4Pe~ 29 ^ 

££ r e =G 7 + 7 S rr -- 7 S] r e(l + l 3)[S r e/(?]+ 2 £ r g. (II.16c) 



III. BOUNDARY CONDITIONS 
A. Liquid crystal analogy 

To motivate the possible boundary conditions in micellar systems, we first recall the situation of nematic liquid 
crystals confined between walls, in which the boundary conditions are a balance between surface and bulk interactions 
[38|. The free energy depends on the nematic order parameter Q = (uu) — ^1, where the unit vector u is parallel to 
the liquid crystalline mesogen [38|. A simple free energy functional for a nematic liquid crystal is [38] 

^lc = / [/ft(Q) + K(VQ) 2 ] dV + \W Q f Tr(Q - Q Q ) 2 dS, (III.l) 

J V J walls 



where /^(Q) is the homogeneous part of the bulk free 
energy, typically expanded in Landau form as /^(Q) = 
|TrQ 2 + . . .; and JC is a Frank elastic constant. In the 
surface free energy, W Q specifies the strength of the wall 
potential which penalizes deviations away from a specific 
order parameter tensor Qo. More complex surface free 
energies are possible (39j , and one may alternatively con- 
sider boundary conditions that influence the orientation, 
but not magnitude, of Q at the wall. 

The spatial dependence of the order parameter is found 
by demanding that the free energy functional be station- 
ary with respect to varying Q. Stationarity in the bulk 
leads to the following differential equation, 



LC 



8Q 



[aQ + ...}- JCV 2 Q = 0, 



(III.2) 



whose boundary condition arises from requiring zero vari- 
ation at the wall: 



JCh- VQ + W Q (Q-Q ) = 0. 



(III.3) 



Here n is the outward unit normal from the fluid at the 
wall. Physically the boundary condition is a balance of 
surface and bulk mesogen torques at the wall. 

The bulk equation defines a correlation length £ = 
yX/* 1 ! which controls the decay of order parameter fluc- 
tuations away from the surface, while the boundary de- 
fines an extrapolation length 



Wa 



(III.4) 



which controls wall anchoring [38(. The extrapolation 
length £ roughly sets the scale for the gradient V ~ l/£ 



induced at the wall by the boundary conditions, in the 
absence of bulk Frank stresses. For small £ the effective 
gradient at the wall is very large, so the wall potential has 
a strong effect. Conversely, for large £ the characteristic 
gradient is small, and the anchoring potential has a weak 
effect. 

The extrapolation length should thus be compared 
with the bulk healing length or correlation length I to 
assess the relevant regime. For strong anchoring (£ <C t) 
the boundary conditions effectively impose Qo at the 
walls, while for weak anchoring (£ 3> £) the boundary 
conditions are dictated by the Frank elastic term, result- 
ing in Neumann boundary conditions, or VQ = 0. Mixed 
boundary conditions apply between these two extremes 

M, 



B. Viscoelastic stress 

Based on the liquid crystalline example, we apply simi- 
lar boundary conditions to the viscoelastic stress; the de- 
tails are given in Appendix [Al The viscoelastic stress is 
analogous to the liquid crystalline order parameter, and 
the stress diffusion constant is analogous to the Frank 
clastic constant (according to T>t — > JC), and there is a 
corresponding anchoring potential, leading to 



Vrn- VS + VF(S-S ) = 0. 



(III.5) 



In this formulation the anchoring potential has the di- 
mensions of length. 

As shown in Appendix the anchoring potential W 
can be expressed as W = Wq/G, where G is the modulus 
and Wo, with dimensions of energy per area, penalizes 
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deviations of the polymer deformation A from a refer- 
ence strain A . A simple, albeit weak, contribution to 
the anchoring potential comes from the effect of the wall 
on nearby micellar or polymer conformations. In Ap- 
pendix [A] we estimate W ~ 2i? s , where R g is the radius 
of gyration. 

By analogy with the liquid crystal example above, the 
extrapolation length is given by 



(III.6) 



IV. METHOD OF CALCULATION 
A. Non-dimensional parameterization 

For both models we express all stresses in units of 
the modulus G, and express time in units of the mi- 
cellar (polymer) stress relaxation time r. In cylindrical 
Couette flow the natural scale for length is the quan- 
tity Riq, which for small q is identical to the gap size 
i?2 — Ri = Ri(l — e~ q )- Hence the dimcnsionless quan- 
tities are: 



where the interfacial width I is set by the diffusion con- 
stant according to I = \JVt . The two characteristic 
lengths are comparable when £ c /£ ~ 1, or as noted 
above, we expect weak anchoring (large £) to yield ef- 
fectively Neumann boundary conditions (zero gradient) 
fl9l . I20I |2(| . Alternatively, very strong anchoring might 
be encouraged by specific wall treatments such as rub- 
bing or grooving the walls along a particular direction 
m, i.e. specifying So ~ mm. This would give very 
small £ and effectively Dirichlet boundary conditions, as 
used by Cook and co-workers [2^, [|(| • 

To determine the extrapolation length an estimate for 
the diffusion constant V is necessary. Two physical ef- 
fects that can lead to non-local dynamics and hence stress 
diffusion are (1) the semi-flexibility of wormlike micelles 
[22| and (2) diffusion of micelles that carry stress [io| . 
In the former case the characteristic length scale is set 
by the persistence length £ p . In Appendix IB"! wc estimate 
these two contributions to be 



, (semi-flexiblility) 
Vt = { 126 v ; 

DtrT (translational diffusion), 



(III.7) 



where _D tr is the translational diffusion coefficient. 

If we consider effects due to semiflexibility, combined 
with the estimate for W due steric wall interactions, we 
estimate anchoring length to be 



\2QR a 



(III. 



For typical giant micelles for which L ~ 200 fim, l v 



670nm, or £ ~ 0.05A. The 
0.094 ^ f 80A, which 



20 nm, we find R g ~ 18/ 
interfacial width in this case is 
would imply strong anchoring. 

Alternatively, if we naively apply the estimate for the 
contribution of translational diffusion to scmidilute mi- 
cellar solutions with D tr ~ 10~ 8 cm 2 /s, r ~ 1 s, one finds 
Vt ~ I/im 2 , leading to £/£ ~ 0.8. Hence in this case 
one might expect the effective boundary conditions to be 
somewhere between Neumann and Dirichlet. 

We stress that these estimates are necessarily very 
crude, and should ultimately be replaced by more pre- 
cise calculations. 



t = t/r 
£ 
G 
W 
Riq' 



£ = 
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(IV. la) 
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(IV.lc) 



In cylindrical Couette flow the boundary conditions are 

. d 

-V— £ + VF(£-£ ) = 0atx = (IV.2) 
ox 

d 

Ve- q — £+W(£-£ ) = 0atx=l, (IV.3) 
ox 

where the change in sign arises from the orientation of 
the surface normal n. 

We use e = 0.05, a = 0.3 for the DJS model and e = 
0.01, j3 = for the DRP model, which places the latter in 
a non-monotonic regime of parameter space (Fig. [J). For 
the shear rate startup calculations shear rates of (7) =3.8 
(DJS model) and (7) =9 (DRP model) were used, which 
are approximately in the middle of the respective stress 
plateaux. 



B. Numerical Methods 



We solve the creeping flow equation (Eq. III.8[) for 
the steady state banding profiles of the two models 
(Eqs. III. 121 and III.16|) for different imposed boundary 
conditions and parameter values. In all cases an aver- 
age shear rate (7) is imposed, so the local shear rate is 
eliminated using Eq. (|II.8|) to obtain a set of coupled sec- 
ond order partial differential equations for E rr , and 
£ r g. This set of equations is then solved subject to the 
chosen boundary conditions by evolving them using the 
Crank-Nicolson algorithm 41 1 . 

Two protocols were followed; the first for shear rate 
sweeps to calculate flow curves, and the second to cal- 
culate shear rate startups from rest. In the first proto- 
col, shear rate sweeps, an initial viscoelastic stress pro- 
file was generated by a set of the 20 longest wavelength 
Fourier modes that fit the desired boundary conditions. 
Randomly assigned amplitudes for modes of the stress 
components S Q( 3 were taken from a uniform distribution 
such that the maximum total value for any component 
1 . The equations of motion were then evolved to find the 
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steady state for the first average imposed shear rate in 
the sweep, which was on either the low or high shear rate 
flow branch of the constitutive curve. This state was then 
used as the initial condition for the next average shear 
rate in the sweep (either incremented or decremented). 
This process was continued until the desired region of 
the flow curve had been calculated. At each shear rate 
the system was evolved typically 500 relaxation times r, 
using a time step of 0.05; it was checked that the steady 
state results did not change for smaller time steps. 

In the second protocol a range of initial viscoclastic 
stress configurations were subjected to a given step in 
the shear rate from zero, and evolved to ensure that the 
same steady state was reached. Using random initial 
conditions often produced multiple bands in this proto- 
col; these multiple bands often did not anneal into two 
bands, which is probably a consequence of both the slow 
motion of interfaces in one dimensional systems and the 
fact that multiple interface solutions are known to be 
locally stable when the interface width is much smaller 
than the system size [42| ■ Because we are not interested, 
here, in the details of this degeneracy, we used a spatially 
smooth initial condition for the viscoelastic stress tensor 
which would encourage the simplest band configuration 
commensurate with the given boundary conditions. A 
variety of different smooth initial conditions were used 
to ensure that the results obtained were robust. 

For smaller V and q steady state could take as long as 
lOOOOr to attain. In the shear banding regime the po- 
sitions of the interfaces between bands were monitored 
to ensure that they had stopped moving. In nearly flat 
geometries with weak curvatures and/or sharp interfaces 
(q < 10~ 4 - 5 ,:D < 10~ 3 ), it was difficult to reliably de- 
termine the steady state interfacial position; this slow 
dynamics of fronts in one dimension is well known [43| . 
Hence we did not study this range of parameter space for 
this protocol. 



The order parameter may be illustrated by plotting an 
ellipsoid with principal axes parallel to and proportional 
to those of the tensor (Fig. [3|. 

By symmetry, we expect the principal axes of So to 
coincide with the axes chosen, although in principle one 
could include non-trivial tilt angles with respect to the 
surface, as is found in the behavior of liquid crystals at 
some interfaces. One could also consider, as was done 
in Refs. [25l.l26j|. a preferred direction of alignment along 
the wall. Below we will consider the natural principal 
axes set by the wall, the principal axes determined by 
the steady state bulk solutions, and a restricted set of 
angles rotated with respect to the natural axes. 
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FIG. 3: Tuning between (top) prolate and oblate boundary 
alignment, and (bottom) homeotropic and tangential align- 
ment. In (top) the principal axes remain unchanged and the 
parameters S and b are varied, while in (bottom) S and b are 
fixed to a prolate distribution, and the angle of orientation 
with respect to the surface is changed. 



C. Notation for (Dirichlet) Boundary Conditions 

A convenient parameterization of the Dirichlet bound- 
ary condition So arises from the physical interpretation 
of viscoelastic stress tensor in terms of the second mo- 
ment of unit vector orientations, S = G(3(uu) — I). In 
its principal frame So can be written (in dimcnsionlcss 
form) as 



iS + |I = i^-b , (IV.4) 
V l -^+b) 

where S = § (cos 2 9) — \ and b = (sin 2 9 cos 2<j>) . The con- 
ventional spherical polar coordinates 9 and <j> describe the 
orientation of unit vectors with respect to the principal 
axes. The parameter S specifics the degree of order along 
the principal direction, and b specifies the degree of bi- 
axiality; S and b obey -\ < S < 1, so <b < 



V. REVIEW: NEUMANN BOUNDARY 
CONDITIONS 

We first review the results for Neumann boundary con- 
ditions, calculated previously for the DJS model [2fJ. In 
this case the diffusion term selects the value T* of the to- 
tal shear stress plateau [24[ . At other values for the total 
stress the interface moves at a speed c ~ T xy — T* y [44| . 
In a flow geometry such as cylindrical Couette flow, in 
which the total stress is inhomogeneous in steady state, 
the interface lies at the position of the flow cell at which 
the total stress is equal to T* y (as long as the stress gra- 
dient is negligible over the scale of the interfacial width 

i = \ff>) [3, [45} . Hence the stress gradient drives the in- 
terface to the correct position. In a flat geometry with no 
stress gradient the speed of approach to the steady state 
position commensurate with the imposed average shear 
rate is determined by interaction with the distant wall, 
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FIG. 4: Hysteresis observed in the DJS model in Couette flow 
for q — 0.00995 upon increasing and decreasing the average 
shear rate (or the gap velocity) along the either edge of the 
stress plateau stress plateau (From [201]). 



and vanishes in the limit of a small diffusion coefficient 

m 

The flow curve in cylindrical Couette flow has a stress 
plateau with a slope dT* /d(j) ~ e q — 1, which vanishes 
in the planar limit q = 0. In either case the flow curves 
display hysteresis at either edge of the stress plateau 
(Fig. [3J [201 . During a shear rate sweep from rest the 
fluid follows the low shear rate branch to a total wall 
stress T r g > T* , until a high shear rate band forms and 
the stress decreases onto the stress plateau. For a down- 
ward shear rate sweep from on the stress plateau the high 
shear rate band shrinks until the interface "touches" the 
wall; at this point a banding solution is unstable with 
respect to a homogeneous phase with a higher stress on 
the low shear rate branch. The width of the hysteresis 
loop decreases upon decreasing the diffusion constant T>. 
Hysteresis of this sort is frequently seen in wormlike mi- 
cellar solutions at the low shear rate edge of the stress 
plateau [l(| HH, E3] ■ The nature of nucleation of a shear 
banded state from a homogeneous phase is unknown and 
worthy of study in its own right. Similar behavior is pre- 
dicted near the high shear rate edge of the stress plateau; 
this has rarely been studied wormlike micelles, because 
in many cases the high shear rate band is unstable or 
unattainable. 

Fig. [5] shows the position of the interface between 
bands function of curvature q, defined as the po- 
sition of midpoint between the high and low values of 
£ r g . The high shear rate band appears to shrink at high 
curvatures (large q) because of the total stress gradient; 
this larger gradient effectively concentrates the shear rate 
near the inner wall, so that the shear band occupies a 
slightly narrower region. The viscoelastic stress tensor 
of the high shear rate band of the DRP model, while 
strongly ordered, is only slightly more aligned with the 
velocity direction than in the low shear rate state. The 
high shear rate state has an alignment angle of 26° which 
is comparable to that reported in experiments on worm- 
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FIG. 5: Position x of the interface between the shear bands 
for Neumann boundary conditions for the DJS (top) and DRP 
(bottom) models, for T> — 1.6 x 10 -3 . The ellipses show the 
eigenvalues and principal axes of the viscoelastic stress tensor, 
which in all cases is nearly prolate in the velocity gradient 
plane. The coordinate x is parallel to the velocity gradient, 



like micelles [lfj. This contrasts with the DJS model, 
in which the viscoelastic stress tensor is strongly aligned 
with the flow in the high shear rate band. Hence the 
DRP model may be a more realistic model for wormlikc 
micelles than the DJS model. 

The upper abscissas of Fig. [5] shows the stress dif- 
ferences in terms of geometric parameters for cylindri- 
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cal Couette flow, AR/R and an equivalent cone angle 
9 for cone and plate flow. In Couette flow, AR/R = 
(i?2 — Ri)/R\ = 1 — e~ q , and the relative stress differ- 
ence between the two cylinders is AT r g/T r g = 1 — e~ 2q . 
In cone and plate geometry the relative stress difference 
is a weak function of the cone angle, ATg^/Tg^, = tan 2 9 
[47l ] . Hence, the cone angle 9 that is roughly equivalent 
to a given Couette curvature satisfies tan 2 9=1 — e~ 2q . 

For reference, the selected stresses for the two models 
in a flat geometry, in the limit of an infinite system and 
Neumann boundary conditions, are: 

—^-= 0.62431 DRP model (V.l) 

f; g y/\ - a 2 = 0.4829 DJS model. (V.2) 



VI. RESULTS 



of this lubricating layer. Hence the effect of the bound- 
ary condition is similar to heterogeneous nucleation or a 
wetting: it provides a site upon which the shear band can 
easily grow, and can eliminate the hysteresis so that the 
shear band grows smoothly from the wall without requir- 
ing "nucleation" . Hysteresis is only seen at the high shear 
rate side of the stress plateau, as shown in profiles 3 and 
6, which are at the same imposed shear rate but after 
different histories. This behavior should be contrasted 
with Neumann boundary conditions, in which hysteresis 
was predicted at both edges of the stress plateau. 

In the stress plateau region the high shear rate band 
lies near the inner cylinder (near x = 0), for both increas- 
ing and decreasing shear rate sweeps (profiles 2 and 5). 
Hence the higher total stress near the wall induces pref- 
erential growth of the inner cylinder's lubricating layer 
upon increasing the average shear rate, while the low 
shear rate band preferentially forms at the outer wall 
upon decreasing the average shear rate. 



We first consider weakly curved Couette flow, and 
study the effect of different Dirichlct boundary conditions 
and the magnitude of the diffusion coefficient, or equiv- 
alently the widths of interfaces, on the flow curves and 
attendant hysteresis. Then we study how the stress pro- 
files and configuration of bands depends on the competi- 
tion between the stress gradient due to the curvature of 
Couettte flow and the boundary conditions. Finally, we 
consider mixed boundary conditions, in which the value 
at the wall can adjust depending on the competition be- 
tween boundary parameters and the stress gradient. 



A. Effect of different Dirichlet boundary conditions 
for fixed geometry 

1. Prolate and flow- aligning anchoring 

Here we choose the degree of curvature to be q = 0.005, 
which corresponds to radii R2 — 1.005i?i, and enforce 
So, the value of the viscoelastic stress at the boundaries 
(Dirichlct boundary conditions). The principal axes of 
So are chosen to be parallel to the flow, gradient and 
vorticity directions. The parametrisation of Eq. (|IV.4|) 
was then applied so that S is completely aligned in the 
flow direction, tangential to the walls, for S — > 1. 

Figure [5] shows the flow curves and shear rate pro- 
files for upward (>) and downward (<) sweeps of the 
DJS model for flow-aligning prolate boundary conditions, 
S = 1.0, b = 0, similar to the value S# of the high shear 
rate branch. For low shear rates the stress is lower than 
that of the low shear rate constitutive curve; profiles I 
and 4. This can be traced to a lubricating layer (evident 
in the inset of profile I) induced by the wall, which has 
a higher shear rate and hence reduces the stress for a 
given imposed shear rate. Hysteresis is not seen at the 
low shear rate edge of the plateau, presumably because 



2. Homeotropic anchoring 

The converse behavior is observed when the bound- 
ary conditions are less favorable to the high shear rate 
phase. Fig. [7] shows the flow profiles and flow curves for 
S = —0.5,6 = —0.5, which corresponds to homeotropic 
orientation with principal axis in the gradient direction. 
In this case a hysteresis loop is found at the low shear 
rate end of the plateau (profiles I and 4), and not at the 
high shear rate end (profiles 3 and 6). This boundary 
condition induces a more viscous layer of low shear rate 
material near the wall, even well into the stress plateau 
(profiles 2 and 5). Profiles 2 and 5 are not quite su- 
perposable, indicating that steady state was probably 
not exactly attained for these profiles; indeed, this value 
T> = 4.4 x 10 -4 is at the lower limit of our computational 
capabilities. The more viscous layer near the boundary 
increases the total stress on the high shear rate branch 
above than that of the constitutive curve (profiles 3 and 
6) , as opposed to the lubricating layer found for the flow- 
aligning boundary condition. 

The shear rate adjacent to the wall is extremely high 
(inset to profile I), even higher than that of the high 
shear rate phase. This is because the principal axes of 
the boundary condition So are aligned with the velocity, 
gradient, and vorticity, such that the shear component 
Sore vanishes. To compensate for this vanishing con- 
tribution to the shear stress, a very high shear rate is 
required to obtain the necessary value for the total shear 
stress. For boundary conditions conducive to the low 
shear rate phase, the viscoelastic stress tensor S adja- 
cent to the wall nonetheless "heals" to the value char- 
acteristic of the low shear rate branch, whereas in the 
flow-aligning case it heals to the value characteristic of 
the high shear rate branch. Although our numerics can 
adequately resolve the sharp gradient near the wall, it 
is likely that higher order gradients are necessary for a 
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FIG. 6: Top left: constitutive curve (solid line) and flow curves (triangles) for upward and downward sweeps of the DJS model 
with T> = 4.4 x 10~ 5 , in cylindrical Couette flow with q — 0.005. The boundary condition is So = diag(2, — 1, —1), as illustrated 
by the ellipsoids (top right). For the noted average imposed shear rates, figures (l)-(6) show the steady state shear rate profiles 
■y(x), and cross-sections of the ellipsoids corresponding to S(x) +1 in the velocity-flow gradient plane. In all cases the ellipsoids 
are nearly prolate. The inset in (1) shows the detail of the shear rate profile near the inner wall x = 0. 



constitutive model to give physically accurate results. 



3. Variation with anchoring angle and from prolate to 
oblate anchoring 

A summary of the behavior of the DJS and DRP mod- 
els for a range of boundary conditions is shown in Fig. [8] 
There are essentially three types of behaviors: hystere- 
sis loops can be observed on (1) the high shear rate end 
of the stress plateau, (2) the low shear rate end of the 
stress plateau, or (3) both ends of the stress plateau. Hys- 
teresis vanishes on those flow branches whose viscoclastic 
stress, loosely, is similar to the imposed boundary value 



So- The two models have broadly similar qualitative be- 
havior. [The apparently greater slope across the plateau 
for the DRP model (Fig. 18 . 2[) is due to the curvature of 
Couette flow and the smaller vertical scale.] 

The behavior was similar when the boundary condi- 
tions were altered by choosing a prolate So and rotating 
its principal axis from the flow gradient (homcotropic) 
to the velocity (tangential) directions. For tangential 
boundary conditions (resembling the high shear rate 
state) there was no overshoot at the start of the plateau 
either on ramping the shear rate up or down. However 
on ramping down the shear rate the flow curve under- 
shoots the stress plateau at the high shear rate end of 
the plateau. As the boundary condition is rotated the 
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FIG. 7: Top left: constitutive curve (solid line) and flow curves (triangles) for upward and downward sweeps of the DJS model 
with T> = 4.4 x 10~ 5 , in cylindrical Couette flow with q = 0.005. The boundary condition is So = diag(— 1, 2, —1), as illustrated 
by the ellipsoid (top right). (l)-(6) show the shear rate profiles 7, and cross-sections of the ellipsoids corresponding to E(x) +1 
in the velocity-flow gradient plane. Except for very near the boundaries, the ellipsoids are nearly prolate. The inset in (1) 
shows the detail of the shear rate profile near the inner wall x = 0. 



hysteresis at the high shear rate end vanishes, and starts 
to develop at the low shear rate end. For homeotropic 
alignment (resembling the low shear rate state) the hys- 
teresis completely vanishes at the high shear rate end of 
the plateau, and hysteresis develops at the low shear rate 
end. 



rate phases is crucial, we focus for the remainder of this 
paper on boundary conditions for which So is one of these 
values. For the DJS model with e = 0.05 and a = 0.3 
these values are 



4. Effect of the diffusion constant 



Having determined that similarity of the imposed vis- 
coelastic stress So to that of either the high or low shear 



12 




iVxv 




-0.22 - 




iogio^) 

.1) Parameter space and characteristic flow curves: DJS Model. 



e) 



0.4 H 
0.2 



-0.2 
-0.4 H 





iogio(i:> 

.2) Parameter space and characteristic flow curves: DRP Model. 



FIG. 8: (a,e) Regions of S — b parameter space for So that give rise to different signatures of hysteresis during shear rate 
ramps, for q — 0.005, and T> = 4.4 x 10~ 5 . The dashed lines on (a,e) correspond to uniaxial order interpolating between an 
oblate "pancake" and a prolate "needle" (at S — —0.5,6 = ±0.5, and S = 1), aligned along the vorticity, flow gradient, and 
flow directions. Hysteresis occurs at the low shear rate end of the stress plateau ( A. b,f); at the high shear rate end of the 
stress plateau (*, d,h); or at both ends shear rate end of the stress plateau (o, c,g). 



13 




FIG. 9: Constitutive curve and flow curves for the DJS model, for q = 0.005, and diffusion constants V = 1.0 x 10 -4 (solid), 
1.6 x 10" 5 (dashed) and 4.0 x 10 -6 (dotted). (a,c) show the start of the plateau, while (b,d) show the end; and the boundary 
condition So equal to the viscoelastic stress Eh or El of, respectively, the high (a,b) or low (c,d) shear rate branch. The 
arrows indicate the direction in which the hysteresis loops are encircled (b,c) or the shear rate swept (a,d). 



S L : (E Lee , S Lrr , S Lr9 ) = (0.40, -0.21, 0.47) low shear rate branch (VLla) 

H H : (Sffee, IW, Sffre) = (1-4, -0.75, 0.12) high shear rate branch, (VI. lb) 

and for the DRP model with (3 = and e — 0.01 they are 

S L : (E Lee , S Lrr , £ Lr9 ) = (0.81, -0.40, 0.60) low shear rate branch (VI.2a) 

T, H : (Snee, SW, Sffre) = (1-7, -0.85, 0.44) high shear rate branch. (VI.2b) 



The extent of the hysteresis depends on the magnitude 
of the diffusion constant. Fig. [9] shows how T>, or equiv- 
alently the width £ ~ spD of the interface, influences 
the flow curves, for the two different boundary condi- 
tions of the DJS model. As the diffusion constant is 
reduced the heterogeneous flow curve approaches the ho- 
mogeneous constitutive curve, and the degree of hystere- 
sis decreases. This occurs for both the DRP (not shown) 
and DJS models. As noted above, this is because the 
interfacial layer occupies a progressively smaller fraction 
of the sample: the separation between the heterogeneous 
flow curve and homogeneous constitutive curve is propor- 
tional to VV ~ L This behavior should be contrasted 
with the Neumann case (Fig.[4]), in which the flow curve 



matches the constitutive curve for all T>. 

B. The effects of flow geometry and total stress 
gradient 

1. Dirichlet boundary conditions 

For a flow cell with a stress gradient, such as a cylin- 
drical Couette cell, the interface in simple banding flow 
with Neumann boundary conditions lies at the position 
in the flow cell where the stress is equal to the selected 
total shear stress T* s . With Dirichlet boundary condi- 
tions there are typically three bands and two interfaces, 
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because of the boundary layer imposed by the wall. We 
now study how the position of these bands varies as a 
function of the flow geometry and diffusion coefficient. 
Calculations were performed using the second protocol 
outlined in Sec. lIVI starting up from rest with a suitable 
initial condition and evolving at fixed average strain rate 
to steady state. The boundary condition So was chosen 
to be either S^ or S#. 

In a flat geometry the phase near the wall is that most 
similar to the imposed boundary conditions, ideally re- 
sulting in a symmetric three band configuration stabi- 
lized by the very weak interactions of the interfaces with 
the walls. For strong enough curvature we expect the 
high shear rate phase to preferentially occupy the high 
stress regions. This is indeed what we find. For increas- 
ing curvature the three band configuration becomes in- 
creasingly asymmetric, as seen in both the DJS and DRP 
models (Figs. HUland [TT|) . A crossover from three banded 
to nearly two-banded behavior occurs at smaller q values 
(weaker curvature) for smaller T>. The three band state 
is most pronounced when the boundaries induce the low 
shear rate phase So- 

The behavior described above can be rationalized as 
follows. For a symmetric flat system (q = 0) both in- 
terfaces lie at the selected stress, since the total stress is 
uniform. The very weak effect a slightly curved system 
q > breaks the symmetry of the shear banded structure, 
and only one (or neither) interface can lie at the selected 
stress; at the same time, the stress gradient would fa- 
vor locating the high shear band near the inner cylinder. 
These two tendencies adjust the band configuration such 
that high (low) shear rate bands grow (shrink) near the 
inner cylinder and shrink (grow) near the outer cylinder. 
For example, for boundary conditions that favor the low 
shear rate band the band near the inner cylinder should 
become narrower with increasing curvature (as seen in 
Fig. [TOa.c). However, a larger diffusion coefficient broad- 
ens the interfaces between bands, which means that the 
walls can strongly influence the bulk behavior. Corre- 
spondingly, for large T> the position of the center shear 
band moves more smoothly as a function of the geometric 
curvature q. For example, compare Fig. [TUb (D = 1CP 2 ) 
with Fig. [TOk (V = 1CT 3 ). 

Just as for Fig. [5J the effective cone angles are again 
shown in the upper abscissas in Figures [TO] and [TTJ For 
small diffusion constant, Fig. [POk shows that in cone and 
plate flow with a typical angle of 9 = 4° , a three banded 
structure could be observed, while the third band would 
be much smaller in Couette flow with AR/R = 0.05, 
and could easily be interpreted as a two band structure. 
Thus, a three band structure, as observed by Britton and 
Callaghan [fj, could be consistent with boundary condi- 
tions that favor the low shear rate form of the viscoelastic 
stress. 

Finally, we note that these results are consistent with 
those of Cook and co-workers [HI, [26[ , who studied a two- 
fluid model of wormlike micelles. In cylindrical Couette 
flow with q = 0.064,1? = 10~ 3 they also found a three- 



band state, but did not explore the effects of geometry 
or diffusion coefficient. 



2. Competing boundary conditions 

Finally we allow surface anchoring to compete with 
the bulk spatial gradients, according to Equation (|III.5[) . 
By varying W, or cquivalently the extrapolation length 
£ = Dt/W, the boundary condition at the wall can be 
tuned smoothly from a fixed value So as W — > oo (£ — > 0) 
to zero gradient as W — ► (£ — > oo). As above, we con- 
sider boundary conditions in which So is the viscoelastic 
stress of either the low or high shear rate branch, S^ or 
S# . We vary both the extrapolation length £ and the ge- 
ometric curvature q, and monitor the behavior of a well- 
formed shear band, with an imposed average shear rate 
such that the high and low shear rate bands are roughly 
the same size. We find qualitatively similar results for 
both the DJS and DRP models. 

Figure [TJ] shows the results for the DRP model for 
boundary conditions So = S^. Consider changing 
the anchoring strength at fixed q (profiles {i,ii,iii,iv} in 
Fig. IT2|) . For weak anchoring strength (large £, profile 
a), the effective boundary condition at the inner wall 
is nearly zero gradient (Neumann), and the high shear 
rate band lies near the inner wall. Upon increasing 
the anchoring strength (reducing £) the boundary con- 
ditions deviates slightly from Neumann conditions be- 
fore reverting to Dirichlet conditions for strong enough 
anchoring (small £), thus imposing S^ at the wall and 
inducing three bands. The crossover occurs when the 
extrapolation length is of order the interfacial thickness, 
£ ~ t = y/Dr. 

In all cases the interface closest to the outer cylinder 
(at larger x) remains at the selected stress T* g while the 
inner interface is at a higher stress, as can be seen by the 
intersections of the shear rate profiles as a function of 
total stress in Fig. irSTiii.iv). Upon approaching the flat 
limit both interfaces approach the selected stress. The 
degree of anchoring strength above which the boundary 
conditions become Dirichlet-likc depends weakly on the 
geometric curvature. For a more highly curved geometry 
a larger anchoring potential (smaller £) is required to en- 
force the boundary condition So = S^. This is because 
the higher stress at the inner wall competes with the ten- 
dency of the wall to induce the low shear rate band. 

We have examined whether or not this "surface tran- 
sition" between different anchoring states has an appre- 
ciable mechanical signal. The torque at the inner cylin- 
der, as would be measured in an experiment, changes 
by of order 1% upon tuning between Neumann-like and 
Dirichlet- like boundary conditions; hence this has only 
a small effect on the macroscopic response of the sys- 
tem. The degree of anchoring is only weakly affected by 
the overall average shear rate; only when one interface 
approaches the wall is there an effect, and again this is 
reflected in only a small change in the measured total 



15 



a) 



AR/R 



0.01 0.03 0.1 0.4 



b) 



i i i i 
e/° 0.4 0.8 1.4 2.5 4.5 8.0 13.9 23.1 34.4 

5j i i i i i i i i i 



AR/R 



0.4 0.8 1.4 2.5 4.5 



0.01 0.03 0.1 0.4 

i 1 1 1 

;.0 13.9 23.1 34.4 




c) 



-5 
AR/R 



-5 
AR/R 



-4 -3 -2 

logio 1 

0.01 0.03 0.1 0.4 

' ' ' 1 di 

0/° 0.4 0.8 1.4 2.5 4.5 8.0 13.9 23.1 34.4 / 

5j i i i i i i i i i 

1 .^ m ^^^^^^^^^^^^^ m ^- 1 



0.8- 



0.6- 



0.4- 



0.2- 




-4 -3 -2 

logio 9 

0.01 0.03 0.1 0.4 

i 1 1 1 

6/° 0.4 0.8 1.4 2.5 4.5 8.0 13.9 23.1 34.4 
i i i i 



0.8- 



0.6- 



0.4- 



0.2 




-3 

logio <i 



FIG. 10: Viscoelastic stress ellipses as a function of position x and curvature q, for the DJS model with (7) = 3.8(log 10 (7) = 
0.58). The diffusion constant is T> — 10 -3 (a,b) and T> = 10 -2 (c,d). The high shear rate phase is lightly shaded (green), and 
the low shear rate phase is dark grey. The boundary conditions are So = Sh or S_l. The top axes show the equivalent Couette 
cylinder gap AR relative to the inner cylinder radius R, and equivalent cone angles 9 if we assume that the stress gradient 
mimics that found in cone and plate flow. 



stress. 



VII. CONCLUSION 



Similar behavior is found for So = S#, with a few 
differences. In this case the interface closest to the in- 
ner wall lies at the selected stress, as can be seen in the 
profiles for j(T r e) in Fig. [13J this implies that for both 
So = Si and So = S# the interface that separates an 
inner band at high shear rate from an outer band at low 
shear rate lies at the selected stress. Another difference 
is that the geometric curvature q has a very weak effect 
on the effective boundary condition, because the higher 
stress at the inner wall doesn't compete with the bound- 
ary condition. 



A. Summary of results 

We have performed a numerical study of the effect of 
the boundary conditions on the viscoelastic (polymer) 
stress for the diffusive Johnson-Segalman (DJS) and dif- 
fusive non-stretching Rolie-Poly (DRP) models, in cylin- 
drical Couette flow. The main results are the following: 

a. Dirichlet boundary conditions (fixed So) influence 
the possible hysteretic behavior for the flow curve: 
for So resembling the high shear rate phase a hys- 
teresis loop occurs at the high shear rate end of 
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FIG. 11: Viscoelastic stress profiles, shown as ellipses as a 
function of position x, for the DRP model in different degrees 
of curvature q, for T> = 10~' ! and (7) = 9 (log 10 (7) = 0.95). 
The high shear rate phase is lightly shaded (green), and the 
low shear rate phase is dark grey. The boundary conditions 
are So = 'Sh (right) or So = Sl (left). The upper abscissa 
shows the equivalent geometrical parameters for cylindrical 
Couette (AR/R) and cone and plate (8) geometries. 



the stress plateau, while for So resembling the low 
shear rate phase a hysteresis loop occurs at the low 
shear rate end of the stress plateau. 

b. The walls can induce a lubrication (or thickened) 
layer that suppresses the hysteresis that would oc- 
cur upon sweeping the shear rate into a banded 
state from a flow branch with characteristic dis- 
similar to the boundary layer. 

c. As with Neumann boundary conditions, the hys- 
teresis loop shrinks with decreasing diffusion con- 



stant. 

d. For Dirichlet boundary conditions the stress distri- 
bution in the gap depends on both the boundary 
condition and the stress gradient. A three-band 
state is stable for weak curvature such as cone and 
plate flow, which becomes more asymmetric or even 
two-banded as the stress gradient increases to that 
of typical cylindrical Couette flow. 

e. For mixed boundary conditions the strength of wall 
anchoring determines whether the effective bound- 
ary conditions are Neumann or Dirichlet. The 
crossover occurs when the extrapolation £ length 
is of order the interfacial width £, or W ~ \J Dr. 

This is the first use of the DRP model to study shear 
banding. The model is microscopically motivated by 
either polymer solutions or wormlike micelles, and its 
non-monotonic behavior arises from better-understood 
physics (convected constraint release competing with 
tube alignment) than that in the DJS model (where the 
"slip parameter" a is necessary). The physics of con- 
vected constraint release yields a less well aligned state 
in the high shear rate phase, whose alignment angle, 
compared to that from the DJS model, is closer to that 
seen experimentally fl6j | . Numerically, interfaces typi- 
cally travel faster in the DRP model than in the DJS 
model, and hence can be more quickly and reliably cal- 
culated, and the DRP model is more robust to large time 
steps and inhomogencous initial conditions. 



B. Discussion and outlook 

An important open problem is the relation of wall slip 
to shear banding flows. A boundary condition such as 
So ~ vVj which imposes alignment parallel to the flow 
direction v, yields velocity profiles that could be inter- 
preted in terms of wall slip, because of the lack of a shear 
component of the viscoelastic stress at the wall. The 
banding profiles in Fig.[6j for example, have small regions 
of very high shear rate near cither wall. Becu et al. mea- 
sured a correlation between band motion and wall slip 
in a solution of CTAB wormlike micelles [H, IH ■ They 
used smooth and sand-blasted Couette cells and found 
different velocity profiles and interface dynamics in the 
banding regime. It is possible that these differences are 
a combination of both different intrinsic boundary con- 
ditions and different degrees of wall slip. 

Another signature of this effect would lie in the mea- 
surement of the low shear rate branch before shear band- 
ing occurs: the flow cell (for example, the smooth one) 
that induced tangential ordering would have a smaller 
stress at a given shear rate, due to the associated lubrica- 
tion layer, than would the flow cell (e.g. the sand-blasted 
one) with less preference for the high shear rate branch. 

In related work, Manneville and co-workers (50j mea- 
sured the rheology and velocity profiles in a triblock 
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FIG. 12: Flow profiles as a function of extrapolation length £ = t>/W and geometrical curvature q for the DRP model with 
mixed boundary conditions, for (7) = 3.8, T> — 7 x 10~ 4 and So equal to the viscoelastic stress in the low shear rate branch, 
Si. (i,iii) Shear rate 7 as a function of position x; the labels a, . . . ,d and 1, . . . , 4 refer to the points in the (£, q) parameter 
map in the center, (ii, iv): Shear rate as a function of the deviation of the total stress T r e ~ e~ 2qx from the selected stress T* . 
(Center) Map of different effective boundary conditions as a function of extrapolation length £ and curvature q: * ~ Neumann 
boundary conditions; A — Dirichlet conditions; o ~ mixed boundary conditions. 



copolymer solution that forms wormlike micelles, using 
heterodyne light scattering. At T = 37° the solutions 
did not shear band but did exhibit wall slip, within the 
~ 40 /im resolution of the technique. At a higher tem- 
perature T = 39.4°, believed to correspond to longer mi- 
celles, the solution exhibited a three-band shear banding 
scenario with the high shear rate bands near the inner 
wall. These data are consistent with the walls inducing 
preferential ordering parallel to the walls with increasing 
temperature and longer micelles, so that that wall slip 
could be either "true" slip or an apparent slip due to the 
specific nature of the boundary conditions. Other mech- 
anisms for slip, such as micelles detaching from the wall 
and disentangling near the surface, may also be relevant 
in this case fBll . [52 . [53| . 

In principle, mixed boundary conditions allows for 
multiple stable inhomogeneous solutions, which could be 
susceptible to non-linear perturbations such as local flow 
inhomogeneities, thermal fluctuations, asperities, motor 
noise, etc. The solutions found here switched between ef- 
fective Neumann and Dirichlet boundary conditions as a 
function of anchoring strength, but easily externally con- 
trollable quantities such as flow geometry and imposed 
average shear rate had very weak effects, and the total 
stress differences between the two kinds of anchoring was 
small (typically less than a percent). 

One could, in principle, envision more complex sur- 
face anchoring behavior, such as a multi-welled wall po- 



tential that would govern a surface phase transition be- 
tween two different wall orientations Soa and £of>- The 
order parameter in the bulk could then, in principle, 
drive the boundary condition between the two poten- 
tial wells, which would then have more dramatic con- 
sequences. Such a bulk-driven wall transition could also 
be triggered by fluctuations, would then lead to erratic 
band motion coupled to an apparent wall slip, depending 
on the details of the specific measurements. 

Another physical effect that has not been considered 
here is the concentration degree of freedom, which in 
some cases probably leads to a depletion layer (or in 
rarer cases perhaps coagulation) layer near the surface. 
Ref. [2(| included concentration in their treatment, but 
without any specific wall-concentration coupling aside 
from a zero flux condition. This is a promising direc- 
tion for future study. 
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FIG. 13: Left: Map of effective boundary conditions as a function of curvature q and extrapolation length £ and selected shear 
rate profiles as a function of position and total stress, for T> — 7 x 10 -4 : * ~ Neumann boundary conditions; A ~ Dirichlet 
conditions; o ~ intermediate boundary conditions. Shear rate as a function of local total stress (middle) and position (right) 
for the DJS and DRP models with boundary conditions So = 'Eh and S = St. 



APPENDIX A: VISCOELASTIC STRESS 
BOUNDARY CONDITIONS 



The dynamics of the DJS or DRP model are formu- 
lated in terms of the viscoelastic (or polymeric) contri- 
bution to the total stress. The free energy of a polymer 
solution is more commonly expressed in terms of the non- 
linear deformation tensor A. A simple free energy func- 
tional that includes surface, bulk, and distortion elastic 



terms is 



F 



d 3 r 



GA 2 + K (VA)' 



d 2 rW (A-A ) 2 
(A.l) 



where K is analogous to the Frank elastic terms that 
penalize distortions of order parameter in nematic liq- 
uid crystals, and G is the modulus. For ordinary poly- 
mer (rubber) elasticity the deformation tensor is A ~ 
(RR) / (Nb 2 ) — hi , where N is the number of Kuhn steps 
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per strand, b is the Kuhn length, and R is the end-to- 
end distance of the strand, and the modulus is thus given 
by G ~ 2>ksTc s , where c s is the number of strands per 
unit volume. A simple dynamics and boundary condi- 
tions for such a free energy (based on the analogy with 
liquid crystals) is 

d t A = -1 [GA - AV 2 A] (A.2) 
= Wo(A-A )+AV 2 A. (A.3) 



Since the viscoelastic stress £ is related to the strain 
by a modulus, S = GA, we can rewrite the equations 
above as 



F = - 
2 



&E = -- 



^ 2 + §(vs) 2 



= Wo (S - £ ) + An • V£, 



(A.4) 

(A.5) 
(A.6) 



where n is the outward unit vector normal to the surface. 
On the other hand, the dynamics of the DJS model is 

1 



[E + X>tV 2 £] 



(A.7) 



from which we identify t = C./G and T>t = K/G. Thus, 
the boundary condition for the DJS model, which follows 
from this correspondence and from Eq. (|A.6j) . is 







X>tV£. 



(A. 



Here Wq is the surface anchoring term that penalizes the 
strain variable at the surface according to Eq. (|A.4|) , In 
the main text we have used W = Wq/G. 



APPENDIX B: ESTIMATES OF W AND V 

One mechanism for such a potential is the steric ex- 
clusion of micelles at the wall, which would favor an 
oblate deformation A and hence an oblate viscoelastic 
stress. Such an effect was found in Monte Carlo sim- 
ulations of polymer melts, which showed a decrease in 
the radius of gyration R^ normal to the wall of 20% 
and a corresponding increase parallel to the wall, inde- 
pendent of the attraction of chain segments to the wall 
[54} . Assuming weak perturbations due to the wall poten- 
tial, the free energy governing subsequent perturbations 
is governed by the entropic elasticity of polymer chains: 
.F/ chain ~ |fceT Tr (A — Aq) 2 . Converting to a wall po- 



tential according to 
Area 



F 

-r-r- c s 2R g 
chain 

3/cbTc s 6 



Tr (A - A ) 



(B.l) 
(B.2) 



we estimate Wq — kBTc s b^/6/N (according to Eq. IA.4|) . 
and W ~ by/2N/3 = 2R g . Here, c s is the number of 
strands per unit volume. 

The stress in a solution of semiflexible polymers, such 
as wormlike micelles, has a contribution S from the tube 
orientation (uu), as well as a liquid crystalline contri- 
bution Slc- Although liquid crystalline effects are only 
appreciable near an isotropic-to-nematic transition, they 
are nonetheless present. Hence, the total stress takes the 
form M,\M 



T = S + 2t?D 



SF lc 



(B.3) 



where F^c is the liquid crystalline contribution to the 
free energy and the neglected terms are non-linear in both 
Q and D. 

For a semiflexible polymer solution Flc has been cal- 
culated by Liu and Fredrickson, and depends on the per- 
sistence length £p. To lowest order in the nematic order 
parameter Q it is 



F, 



LC 



[aTr(Q 2 )+/Ca M Q:9 M Q], (B.4) 



where only a single Frank constant has been included 
here. This leads to a liquid crystalline contribution to 
the stress of the form 



Slc = aQ -/CV 2 Q. 
In this case the parameters are [22[ 

45 k B T£ p (f> 45# B T 



K 



126tt D 2 



nD 2 £ 



(B.5) 



(B.6) 



where <p is the volume fraction <p [221 ] and D is the micellar 
diameter. For wormlike micelles £ p ~ 20 nm, D ~ 2 nm 
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[HI, leading to JC ~ 0.37 <fi ieT/nm. Upon converting 
the entire description to one in terms of the total micel- 
lar stress S tot = S + S^c" an d performing a gradient 
expansion, one finds the following contribution to the 
stress diffusion coefficient due to Frank-like elasticity: 



There may of course, be other non-local contributions, 
which could depend on micellar size, concentration, hy- 
drodynamics, or other physics. 
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